# This file contains the S-plus code that was used to implement the estimators and functionals in "Histogram for Hazard Rate Estimation", Sankhya B, (74): 286-301, (2013).
# As a fully working example the code herein can reproduce Figure 1b of the paper (the relevant data are also given below). 
# The same functions (with different input data) can be used in reproducing other figures of the same paper 
# More general use of the estimates is also feasible by supplying the corresponding data set.  Use of the functions is self - explanatory and is further facilitated by the comments at he preamble of each function.

#N.B. Compatibility with R is not guaranteed as the functions were not written for use with R - however, with minimal modifications, their use in the R environment should be feasible.


"lognormder"<-function(x, sigma)
{
	( -1/(x^2 *  sigma * sqrt(2*pi) )    * exp( - ((logb(x, base=exp(1)))^2 )/(2*sigma^2) ) * (1- logb(x, base=exp(1))/sigma^2  ))^2
}

"SimpsonInt"<- function(xin, h)
{   #xin = data, h=distance of the points at which function is evaluated
    n<-length(xin)
    nn<-1:n
    int0<-xin[1]
    int2n<-xin[n]
    inteven<-xin[seq(3,n-1,by=2)]
    intodd<-xin[seq(2,n-2,by=2)]
    sumodd<-sum(intodd)
    sumeven<-sum(inteven)
    res<-(h/3)*(int0+4*sumodd+2*sumeven+int2n)
    res
}


"lnormhaz"<-function(x, sigma)
{
	dlnorm(x,0,sigma)/(1-plnorm(x,0,sigma))
}


"lognhazder"<-function(x, sig, a)
{
	#a<-lognormder(x, sig)
	((a * (1-plnorm(x, 0, sig)) + dlnorm(x,0,sig))/((1-plnorm(x,0,sig))^2))^2
}


hnstar<-function(x, eta, samplesize)
{
	weibdersq<-lognormder(x, eta)
		h<-x[2]-x[1]
	a<-SimpsonInt(weibdersq, h)
	(6/(samplesize * a))^(1/3)
}


hn0<-function(x, eta, samplesize)
{
	lnormhr<-lnormhaz(x, eta)
	b<-plnorm(x,0, eta)
	h<-x[2]-x[1]
	numintegrand<-lnormhr/b
	num<-SimpsonInt(numintegrand, h)
	dd<-lognormder(x,eta)
	
	d<- lognhazder(x,eta, dd) + 3* lnormhr^4
	#d<-  3* lnormhr^4

	den<-samplesize * SimpsonInt(d, h)
	(6*num / den)^(1/3)
	
}


xout<-seq(0.01, .5, length=100)
etaseq<-seq(0, 2, length = 100)
samplesize<-400

nn<-length(etaseq)
out<-sapply(1:nn, function(i, xout, samplesize, etaseq)
						hnstar(xout, etaseq[i], samplesize)
						, xout, samplesize, etaseq)
						
out1<-sapply(1:nn, function(i, xout, samplesize, etaseq)
						hn0(xout, etaseq[i], samplesize)
						, xout, samplesize, etaseq)

						
plot(etaseq, out1, type="l", xlab="sigma parameter", ylab="optimal bindwidth", ylim=c(0,0.5))
lines(etaseq, out, lty=4)
key(text=list(c("h_n^0", " ",  "h_n^\star"), font =4), 
lines = list(lty=c(1,0,4)), corner=c(1,1),
cex= 1.5)